/*
 * Copyright 1993-2007 NVIDIA Corporation.  All rights reserved.
 *
 * NOTICE TO USER:
 *
 * This source code is subject to NVIDIA ownership rights under U.S. and
 * international Copyright laws.  Users and possessors of this source code
 * are hereby granted a nonexclusive, royalty-free license to use this code
 * in individual and commercial software.
 *
 * NVIDIA MAKES NO REPRESENTATION ABOUT THE SUITABILITY OF THIS SOURCE
 * CODE FOR ANY PURPOSE.  IT IS PROVIDED "AS IS" WITHOUT EXPRESS OR
 * IMPLIED WARRANTY OF ANY KIND.  NVIDIA DISCLAIMS ALL WARRANTIES WITH
 * REGARD TO THIS SOURCE CODE, INCLUDING ALL IMPLIED WARRANTIES OF
 * MERCHANTABILITY, NONINFRINGEMENT, AND FITNESS FOR A PARTICULAR PURPOSE.
 * IN NO EVENT SHALL NVIDIA BE LIABLE FOR ANY SPECIAL, INDIRECT, INCIDENTAL,
 * OR CONSEQUENTIAL DAMAGES, OR ANY DAMAGES WHATSOEVER RESULTING FROM LOSS
 * OF USE, DATA OR PROFITS,  WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE
 * OR OTHER TORTIOUS ACTION,  ARISING OUT OF OR IN CONNECTION WITH THE USE
 * OR PERFORMANCE OF THIS SOURCE CODE.
 *
 * U.S. Government End Users.   This source code is a "commercial item" as
 * that term is defined at  48 C.F.R. 2.101 (OCT 1995), consisting  of
 * "commercial computer  software"  and "commercial computer software
 * documentation" as such terms are  used in 48 C.F.R. 12.212 (SEPT 1995)
 * and is provided to the U.S. Government only as a commercial end item.
 * Consistent with 48 C.F.R.12.212 and 48 C.F.R. 227.7202-1 through
 * 227.7202-4 (JUNE 1995), all U.S. Government End Users acquire the
 * source code with only those rights set forth herein.
 *
 * Any use of this source code in individual and commercial software must
 * include, in the user documentation and internal comments to the code,
 * the above Disclaimer and U.S. Government End Users Notice.
 */

#include <stdio.h>
#include <stdlib.h>
#include <cutil_inline.h>

#include "SobelFilter_kernels.h"

#define SV 0.003921f
#define IV 255.f

//for SIND
#define TSS 2.0f

#define TriMatWidth 512

//__global__  float SNR;
#define BlockSize 16
static float *SNR = NULL;

// Texture reference for reading image
texture<unsigned char, 2> tex;
texture<float, 2> ftex;

extern __shared__ unsigned char LocalBlock[]; 
//for SIND
//extern __shared__ float localROI[];
//extern __shared__ float sqrlocalROI[];

static cudaArray *array = NULL;
static cudaArray *farray = NULL;

int imin;
int imax;

#define Radius 1

#ifdef FIXED_BLOCKWIDTH
#define BlockWidth 80
#define SharedPitch 384
#endif



__device__ unsigned char
ComputeSobel(unsigned char ul, // upper left
             unsigned char um, // upper middle
             unsigned char ur, // upper right
             unsigned char ml, // middle left
             unsigned char mm, // middle (unused)
             unsigned char mr, // middle right
             unsigned char ll, // lower left
             unsigned char lm, // lower middle
             unsigned char lr, // lower right
             float fScale )
{
    short Horz = ur + 2*mr + lr - ul - 2*ml - ll;
    short Vert = ul + 2*um + ur - ll - 2*lm - lr;
    short Sum = (short) (fScale*(abs(Horz)+abs(Vert)));
    if ( Sum < 0 ) return 0; else if ( Sum > 0xff ) return 0xff;
    return (unsigned char) Sum;
}

__device__ float
SIND_Derivative_App(float4 f4Valueneighbor, float fCurvalue)
{
	float fTemp = __powf((f4Valueneighbor.y-fCurvalue), 2) + __powf((f4Valueneighbor.w-fCurvalue), 2) +
	/*float fTempL =*/ __powf((fCurvalue-f4Valueneighbor.x), 2) + __powf((fCurvalue-f4Valueneighbor.z), 2);

	return fTemp;

}

__device__ float
SIND_Laplacian_App(float4 f4Valueneighbor, float fCurvalue)
{
	return ((f4Valueneighbor.x + f4Valueneighbor.y + f4Valueneighbor.z + f4Valueneighbor.w - 4.0f*fCurvalue)/fCurvalue);
}

__device__ float
SIND_Q_Function(float fSqrDerivative, float fLaplacian, float fCurvalue)
{
	float numerator = 0.5f*fSqrDerivative/(fCurvalue*fCurvalue) - 1.0f/16.0f*(fLaplacian*fLaplacian);
	float denominator = 1 + 0.25f*fLaplacian;
	return __fsqrt_rn(float(numerator/(denominator*denominator)));
}


__device__ float
SIND_C_Function(float fSNR, float4 f4Valueneighbor, float fCurvalue)
{
	float fSqrDerivative = SIND_Derivative_App(f4Valueneighbor, fCurvalue);
	float fLaplacian = SIND_Laplacian_App(f4Valueneighbor, fCurvalue);

	float q = SIND_Q_Function(fSqrDerivative, fLaplacian, fCurvalue);

	float output = 1.0/(1.0+(q*q-fSNR*fSNR)/(fSNR*fSNR*(1+fSNR*fSNR)));
	
	return output;
}


__global__ void 
SobelShared( uchar4 *pSobelOriginal, unsigned short SobelPitch, 
#ifndef FIXED_BLOCKWIDTH
             short BlockWidth, short SharedPitch,
#endif
             short w, short h, float fScale )
{ 
    short u = 4*blockIdx.x*BlockWidth;
    short v = blockIdx.y*blockDim.y + threadIdx.y;
    short ib;

    int SharedIdx = threadIdx.y * SharedPitch;

    for ( ib = threadIdx.x; ib < BlockWidth+2*Radius; ib += blockDim.x ) {
        LocalBlock[SharedIdx+4*ib+0] = tex2D( tex, 
            (float) (u+4*ib-Radius+0), (float) (v-Radius) );
        LocalBlock[SharedIdx+4*ib+1] = tex2D( tex, 
            (float) (u+4*ib-Radius+1), (float) (v-Radius) );
        LocalBlock[SharedIdx+4*ib+2] = tex2D( tex, 
            (float) (u+4*ib-Radius+2), (float) (v-Radius) );
        LocalBlock[SharedIdx+4*ib+3] = tex2D( tex, 
            (float) (u+4*ib-Radius+3), (float) (v-Radius) );
    }
    if ( threadIdx.y < Radius*2 ) {
        //
        // copy trailing Radius*2 rows of pixels into shared
        //
        SharedIdx = (blockDim.y+threadIdx.y) * SharedPitch;
        for ( ib = threadIdx.x; ib < BlockWidth+2*Radius; ib += blockDim.x ) {
            LocalBlock[SharedIdx+4*ib+0] = tex2D( tex, 
                (float) (u+4*ib-Radius+0), (float) (v+blockDim.y-Radius) );
            LocalBlock[SharedIdx+4*ib+1] = tex2D( tex, 
                (float) (u+4*ib-Radius+1), (float) (v+blockDim.y-Radius) );
            LocalBlock[SharedIdx+4*ib+2] = tex2D( tex, 
                (float) (u+4*ib-Radius+2), (float) (v+blockDim.y-Radius) );
            LocalBlock[SharedIdx+4*ib+3] = tex2D( tex, 
                (float) (u+4*ib-Radius+3), (float) (v+blockDim.y-Radius) );
        }
    }

    __syncthreads();

    u >>= 2;    // index as uchar4 from here
    uchar4 *pSobel = (uchar4 *) (((char *) pSobelOriginal)+v*SobelPitch);
    SharedIdx = threadIdx.y * SharedPitch;

    for ( ib = threadIdx.x; ib < BlockWidth; ib += blockDim.x ) {

        unsigned char pix00 = LocalBlock[SharedIdx+4*ib+0*SharedPitch+0];
        unsigned char pix01 = LocalBlock[SharedIdx+4*ib+0*SharedPitch+1];
        unsigned char pix02 = LocalBlock[SharedIdx+4*ib+0*SharedPitch+2];
        unsigned char pix10 = LocalBlock[SharedIdx+4*ib+1*SharedPitch+0];
        unsigned char pix11 = LocalBlock[SharedIdx+4*ib+1*SharedPitch+1];
        unsigned char pix12 = LocalBlock[SharedIdx+4*ib+1*SharedPitch+2];
        unsigned char pix20 = LocalBlock[SharedIdx+4*ib+2*SharedPitch+0];
        unsigned char pix21 = LocalBlock[SharedIdx+4*ib+2*SharedPitch+1];
        unsigned char pix22 = LocalBlock[SharedIdx+4*ib+2*SharedPitch+2];

        uchar4 out;

        out.x = ComputeSobel(pix00, pix01, pix02, 
                             pix10, pix11, pix12, 
                             pix20, pix21, pix22, fScale );

        pix00 = LocalBlock[SharedIdx+4*ib+0*SharedPitch+3];
        pix10 = LocalBlock[SharedIdx+4*ib+1*SharedPitch+3];
        pix20 = LocalBlock[SharedIdx+4*ib+2*SharedPitch+3];
        out.y = ComputeSobel(pix01, pix02, pix00, 
                             pix11, pix12, pix10, 
                             pix21, pix22, pix20, fScale );

        pix01 = LocalBlock[SharedIdx+4*ib+0*SharedPitch+4];
        pix11 = LocalBlock[SharedIdx+4*ib+1*SharedPitch+4];
        pix21 = LocalBlock[SharedIdx+4*ib+2*SharedPitch+4];
        out.z = ComputeSobel( pix02, pix00, pix01, 
                              pix12, pix10, pix11, 
                              pix22, pix20, pix21, fScale );

        pix02 = LocalBlock[SharedIdx+4*ib+0*SharedPitch+5];
        pix12 = LocalBlock[SharedIdx+4*ib+1*SharedPitch+5];
        pix22 = LocalBlock[SharedIdx+4*ib+2*SharedPitch+5];
        out.w = ComputeSobel( pix00, pix01, pix02, 
                              pix10, pix11, pix12, 
                              pix20, pix21, pix22, fScale );
        if ( u+ib < w/4 && v < h ) {
            pSobel[u+ib] = out;
        }
    }

    __syncthreads();
}

__global__ void
SIND_SNR(unsigned int Pitch, int w, int h, int r1, int c1, int r2, int c2, float *SNR)
{
	//Thread index
	int tx = threadIdx.x;
	int ty = threadIdx.y;
	int tid = ty*blockDim.x + tx;

	__shared__ float mean;
	__shared__ float variance;
	//__shared__ float sqrsum;
	//__shared__ float localROI[BlockSize][BlockSize];
	__shared__ float localROI[BlockSize*BlockSize];
	__shared__ float sqrlocalROI[BlockSize*BlockSize];
	int threadno = BlockSize*BlockSize;

	//localROI[ty][tx] = (float)pImage[(r1+ty)*Pitch+c1+tx];
	
	float temp = (float)tex2D( ftex, (float) c1+tx, (float) r1+ty );
	//pSobel[i] = min( max((tex2D( tex, (float) i, (float) blockIdx.x ) * fscale), 0.f), 255.f);
	localROI[tid] = temp;
	sqrlocalROI[tid] = temp*temp;

	__syncthreads();
	
	if(threadno >= 1024)
	{
		if(tid < 512)
		{
			localROI[tid] += localROI[tid + 512];
			sqrlocalROI[tid] += sqrlocalROI[tid+512];
		}
		__syncthreads();

	}
	

	if(threadno >= 512)
	{
		if(tid < 256)
		{
			localROI[tid] += localROI[tid + 256];
			sqrlocalROI[tid] += sqrlocalROI[tid + 256];
		}
		__syncthreads();
	}
	
	if(threadno >= 256)
	{
		if(tid < 128)
		{
			localROI[tid] += localROI[tid + 128];
			sqrlocalROI[tid] += sqrlocalROI[tid+128];
		}
		__syncthreads();

	}
	
	if(threadno >= 128)
	{
		if(tid < 64)
		{
			localROI[tid] += localROI[tid + 64];
			sqrlocalROI[tid] += sqrlocalROI[tid + 64];
		}
		__syncthreads();
	}

	
	if(tid < 32)
	{
		if(threadno >= 64)
		{
			localROI[tid] += localROI[tid + 32];
			sqrlocalROI[tid] += sqrlocalROI[tid + 32];
		}
		if(threadno >= 32)
		{
			localROI[tid] += localROI[tid + 16];
			sqrlocalROI[tid] += sqrlocalROI[tid + 16];	
		}
		if(threadno >= 16)
		{
			localROI[tid] += localROI[tid + 8];
			sqrlocalROI[tid] += sqrlocalROI[tid + 8];
		}
		if(threadno >= 8)
		{
			localROI[tid] += localROI[tid + 4];
			sqrlocalROI[tid] += sqrlocalROI[tid + 4];
		}
		if(threadno >= 4)
		{
			localROI[tid] += localROI[tid + 2];
			sqrlocalROI[tid] += sqrlocalROI[tid + 2];
		}

		if(threadno >= 2)
		{
			localROI[tid] += localROI[tid + 1];
			sqrlocalROI[tid] += sqrlocalROI[tid + 1];
		}
	}

	if(tid == 0) 
	{
		mean = localROI[0]/(float)(BlockSize*BlockSize);
		variance = sqrlocalROI[0]/(float)(BlockSize*BlockSize) - mean*mean;
		
		*SNR = __fsqrt_rn(variance)/mean;
		//*SNR = localROI[0];
	}

}

__global__ void 
SobelCopyImage( Pixel *pSobelOriginal, unsigned int Pitch, 
                int w, int h, float fscale )
{ 
    unsigned char *pSobel = 
      (unsigned char *) (((char *) pSobelOriginal)+blockIdx.x*Pitch);
    for ( int i = threadIdx.x; i < w; i += blockDim.x ) {
        pSobel[i] = min( max((tex2D( tex, (float) i, (float) blockIdx.x ) * fscale), 0.f), 255.f);
    }
}

__global__ void
LogCompress(float *pLogCompressedImg, unsigned int Pitch,
				int w, int h, int *mn, int *mx)
{
	float *pLogCompress = 
		(float *)(((float *)pLogCompressedImg) + blockIdx.x*Pitch);
	for ( int i = threadIdx.x; i < w; i += blockDim.x ) 
	{
		//float temp = (float)min( max((float)tex2D( tex, (float) tx, (float) ty ), 0.f), 255.f);
        float temp =(float)min( max((float)tex2D( tex, (float) i, (float) blockIdx.x ), 0.f), 255.f);
		pLogCompress[i] = __expf((temp-*mn)/float(*mx - *mn));

    }
}

__global__ void 
Diffusivity(float4 *pDiffMat, int Pitch, 
			int w, int h, float fSNR)
{
	int tix = blockIdx.x*blockDim.x + threadIdx.x;
	int tiy = blockIdx.y*blockDim.y + threadIdx.y;

	float4 value;
	float4 valueNext;
	float4 valueN;
	float4 valueS;
	float4 valueW;
	float4 valueE;

	float curvalue = (float)tex2D(ftex, (float)tix, (float)tiy);

	value.x = (float)tex2D(ftex, (float)tix, float(tiy-1));
	value.y = (float)tex2D(ftex, (float)tix, float(tiy+1));
	value.z = (float)tex2D(ftex, float(tix-1), (float)tiy);
	value.w = (float)tex2D(ftex, float(tix+1), (float)tiy);

	valueNext.x = (float)tex2D(ftex, (float)(tix+1), float(tiy));
	valueNext.y = (float)tex2D(ftex, (float)(tix+1), float(tiy+2));
	valueNext.z = (float)tex2D(ftex, float(tix), (float)(tiy+1));
	valueNext.w = (float)tex2D(ftex, float(tix+2), (float)(tiy+1));

	valueN.x = (float)tex2D(ftex, (float)tix, float(tiy-2));
	valueN.y = (float)tex2D(ftex, (float)tix, float(tiy));
	valueN.z = (float)tex2D(ftex, float(tix-1), (float)(tiy-1));
	valueN.w = (float)tex2D(ftex, float(tix+1), (float)(tiy-1));

	valueS.x = (float)tex2D(ftex, (float)tix, float(tiy));
	valueS.y = (float)tex2D(ftex, (float)tix, float(tiy+2));
	valueS.z = (float)tex2D(ftex, float(tix-1), (float)(tiy+1));
	valueS.w = (float)tex2D(ftex, float(tix+1), (float)(tiy+1));

	valueW.x = (float)tex2D(ftex, (float)(tix-1), float(tiy-1));
	valueW.y = (float)tex2D(ftex, (float)(tix-1), float(tiy+1));
	valueW.z = (float)tex2D(ftex, float(tix-2), (float)tiy);
	valueW.w = (float)tex2D(ftex, float(tix), (float)tiy);

	valueE.x = (float)tex2D(ftex, (float)(tix+1), float(tiy-1));
	valueE.y = (float)tex2D(ftex, (float)(tix+1), float(tiy+1));
	valueE.z = (float)tex2D(ftex, float(tix), (float)tiy);
	valueE.w = (float)tex2D(ftex, float(tix+2), (float)tiy);

	float curC = SIND_C_Function(fSNR, value, curvalue);
	float nextC = SIND_C_Function(fSNR, valueNext, curvalue);
	float nC = SIND_C_Function(fSNR, valueN, curvalue);
	float sC = SIND_C_Function(fSNR, valueS, curvalue);
	float eC = SIND_C_Function(fSNR, valueE, curvalue);
	float wC = SIND_C_Function(fSNR, valueW, curvalue);

	pDiffMat[tix*w+tiy].x = 0.5*(nC + curC);
	pDiffMat[tix*w+tiy].y = 0.5*(sC + curC);
	pDiffMat[tix*w+tiy].z = 0.5*(eC + curC);
	pDiffMat[tix*w+tiy].w = 0.5*(wC + curC);

}

__global__ void
CyclicReduction(float4 *pDiffMat, Pixel *pOutput, float *pTriImgbk, int Pitch, int w, int h, float lambda )
{
	int tix = blockIdx.x*blockDim.x + threadIdx.x;
	int tiy = blockIdx.y*blockDim.y + threadIdx.y;

	__shared__ float3 TriMat[512];
	__shared__ float3 TriMat2[255];
	__shared__ float3 TriMat3[127];
	__shared__ float3 TriMat4[63];
	__shared__ float3 TriMat5[31];
	__shared__ float3 TriMat6[15];
	__shared__ float3 TriMat7[7];
	__shared__ float3 TriMat8[3];
	__shared__ float3 TriMat9[1];
	//__shared__ float3 TriMat10[1];


	__shared__ float TriImg[512];
	__shared__ float TriImg2[255];
	__shared__ float TriImg3[127];
	__shared__ float TriImg4[63];
	__shared__ float TriImg5[31];
	__shared__ float TriImg6[15];
	__shared__ float TriImg7[7];
	__shared__ float TriImg8[3];
	__shared__ float TriImg9[1];
	//__shared__ float3 TriMat10[1];

	  float *TriImgbk = pTriImgbk+512;
	  float *TriImgbk2 = pTriImgbk+512+255;
	  float *TriImgbk3 = pTriImgbk+512+255+127;
	  float *TriImgbk4 = pTriImgbk+512+255+127+63;
	  float *TriImgbk5 = pTriImgbk+512+255+127+63+31;
	  float *TriImgbk6 = pTriImgbk+512+255+127+63+31+15;
	  float *TriImgbk7 = pTriImgbk+512+255+127+63+31+15+7;
	  float *TriImgbk8 = pTriImgbk+512+255+127+63+31+15+7+3;
	  float *TriImgbk9 = pTriImgbk+512+255+127+63+31+15+7+3+1;

// 	__shared__ float TriImgbk[512];
// 	__shared__ float TriImgbk2[255];
// 	__shared__ float TriImgbk3[127];
// 	__shared__ float TriImgbk4[63];
// 	__shared__ float TriImgbk5[31];
// 	__shared__ float TriImgbk6[15];
// 	__shared__ float TriImgbk7[7];
// 	__shared__ float TriImgbk8[3];
// 	__shared__ float TriImgbk9[1];


	//Row
	int tidx = tiy*w+tix;
	TriMat[threadIdx.x].x = -4.0f*lambda*pDiffMat[tidx].x;
	TriMat[threadIdx.x].y = 2.0f - 4.0f*lambda*(-(pDiffMat[tidx].y+pDiffMat[tidx].x));
	TriMat[threadIdx.x].z = -4.0f*lambda*pDiffMat[tidx].y;

	TriImg[threadIdx.x] = (float)tex2D(ftex, (float)tix, (float)tiy);

	__syncthreads();

	//Solve Tridiagnal System
	float aa, rr;
	//Forward
	if(threadIdx.x < 255)
	{
		aa = TriMat[2*threadIdx.x+1].x/TriMat[2*threadIdx.x].y;
		rr = TriMat[2*threadIdx.x+1].z/TriMat[2*threadIdx.x+2].y;

		TriMat2[threadIdx.x].x = -aa*TriMat[2*threadIdx.x].x;
		TriMat2[threadIdx.x].y = TriMat[2*threadIdx.x+1].y - aa*TriMat[2*threadIdx.x].z - rr*TriMat[2*threadIdx.x+2].x;
		TriMat2[threadIdx.x].z = -rr*TriMat[2*threadIdx.x+2].z;

		TriImg2[threadIdx.x] = TriImg[2*threadIdx.x+1] - aa*TriImg[2*threadIdx.x] - rr*TriImg[2*threadIdx.x+2];
	}
	__syncthreads();

	
	if(threadIdx.x < 127)
	{
		aa = TriMat2[2*threadIdx.x+1].x/TriMat2[2*threadIdx.x].y;
		rr = TriMat2[2*threadIdx.x+1].z/TriMat2[2*threadIdx.x+2].y;

		TriMat3[threadIdx.x].x = -aa*TriMat2[2*threadIdx.x].x;
		TriMat3[threadIdx.x].y = TriMat2[2*threadIdx.x+1].y - aa*TriMat2[2*threadIdx.x].z - rr*TriMat2[2*threadIdx.x+2].x;
		TriMat3[threadIdx.x].z = -rr*TriMat2[2*threadIdx.x+2].z;

		TriImg3[threadIdx.x] = TriImg2[2*threadIdx.x+1] - aa*TriImg2[2*threadIdx.x] - rr*TriImg2[2*threadIdx.x+2];
	}
	__syncthreads();

	if(threadIdx.x <63)
	{
		aa = TriMat3[2*threadIdx.x+1].x/TriMat3[2*threadIdx.x].y;
		rr = TriMat3[2*threadIdx.x+1].z/TriMat3[2*threadIdx.x+2].y;

		TriMat4[threadIdx.x].x = -aa*TriMat3[2*threadIdx.x].x;
		TriMat4[threadIdx.x].y = TriMat3[2*threadIdx.x+1].y - aa*TriMat3[2*threadIdx.x].z - rr*TriMat3[2*threadIdx.x+2].x;
		TriMat4[threadIdx.x].z = -rr*TriMat3[2*threadIdx.x+2].z;

		TriImg4[threadIdx.x] = TriImg3[2*threadIdx.x+1] - aa*TriImg3[2*threadIdx.x] - rr*TriImg2[2*threadIdx.x+2];
	}
	__syncthreads();

		if(threadIdx.x < 31)
	{
		aa = TriMat4[2*threadIdx.x+1].x/TriMat4[2*threadIdx.x].y;
		rr = TriMat4[2*threadIdx.x+1].z/TriMat4[2*threadIdx.x+2].y;

		TriMat5[threadIdx.x].x = -aa*TriMat4[2*threadIdx.x].x;
		TriMat5[threadIdx.x].y = TriMat4[2*threadIdx.x+1].y - aa*TriMat4[2*threadIdx.x].z - rr*TriMat4[2*threadIdx.x+2].x;
		TriMat5[threadIdx.x].z = -rr*TriMat4[2*threadIdx.x+2].z;

		TriImg5[threadIdx.x] = TriImg4[2*threadIdx.x+1] - aa*TriImg4[2*threadIdx.x] - rr*TriImg4[2*threadIdx.x+2];
	}
	__syncthreads();
	
	if(threadIdx.x < 15)
	{
		aa = TriMat5[2*threadIdx.x+1].x/TriMat5[2*threadIdx.x].y;
		rr = TriMat5[2*threadIdx.x+1].z/TriMat5[2*threadIdx.x+2].y;

		TriMat6[threadIdx.x].x = -aa*TriMat5[2*threadIdx.x].x;
		TriMat6[threadIdx.x].y = TriMat5[2*threadIdx.x+1].y - aa*TriMat5[2*threadIdx.x].z - rr*TriMat5[2*threadIdx.x+2].x;
		TriMat6[threadIdx.x].z = -rr*TriMat5[2*threadIdx.x+2].z;

		TriImg6[threadIdx.x] = TriImg5[2*threadIdx.x+1] - aa*TriImg5[2*threadIdx.x] - rr*TriImg5[2*threadIdx.x+2];
	}
	__syncthreads();

	if(threadIdx.x < 7)
	{
		aa = TriMat6[2*threadIdx.x+1].x/TriMat6[2*threadIdx.x].y;
		rr = TriMat6[2*threadIdx.x+1].z/TriMat6[2*threadIdx.x+2].y;

		TriMat7[threadIdx.x].x = -aa*TriMat6[2*threadIdx.x].x;
		TriMat7[threadIdx.x].y = TriMat6[2*threadIdx.x+1].y - aa*TriMat6[2*threadIdx.x].z - rr*TriMat6[2*threadIdx.x+2].x;
		TriMat7[threadIdx.x].z = -rr*TriMat6[2*threadIdx.x+2].z;

		TriImg7[threadIdx.x] = TriImg6[2*threadIdx.x+1] - aa*TriImg6[2*threadIdx.x] - rr*TriImg6[2*threadIdx.x+2];
	}
	__syncthreads();
	
	if(threadIdx.x < 3)
	{
		aa = TriMat7[2*threadIdx.x+1].x/TriMat7[2*threadIdx.x].y;
		rr = TriMat7[2*threadIdx.x+1].z/TriMat7[2*threadIdx.x+2].y;

		TriMat8[threadIdx.x].x = -aa*TriMat7[2*threadIdx.x].x;
		TriMat8[threadIdx.x].y = TriMat7[2*threadIdx.x+1].y - aa*TriMat7[2*threadIdx.x].z - rr*TriMat7[2*threadIdx.x+2].x;
		TriMat8[threadIdx.x].z = -rr*TriMat7[2*threadIdx.x+2].z;

		TriImg8[threadIdx.x] = TriImg7[2*threadIdx.x+1] - aa*TriImg7[2*threadIdx.x] - rr*TriImg7[2*threadIdx.x+2];
	}
	__syncthreads();

	
	if(threadIdx.x < 1)
	{
		aa = TriMat8[2*threadIdx.x+1].x/TriMat8[2*threadIdx.x].y;
		rr = TriMat8[2*threadIdx.x+1].z/TriMat8[2*threadIdx.x+2].y;

		TriMat9[threadIdx.x].x = -aa*TriMat8[2*threadIdx.x].x;
		TriMat9[threadIdx.x].y = TriMat8[2*threadIdx.x+1].y - aa*TriMat8[2*threadIdx.x].z - rr*TriMat8[2*threadIdx.x+2].x;
		TriMat9[threadIdx.x].z = -rr*TriMat8[2*threadIdx.x+2].z;

		TriImg9[threadIdx.x] = TriImg8[2*threadIdx.x+1] - aa*TriImg8[2*threadIdx.x] - rr*TriImg8[2*threadIdx.x+2];

		TriImg9[threadIdx.x] = TriImg9[threadIdx.x]/TriMat9[threadIdx.x].y;
	}
	__syncthreads();

	//Back
	if(threadIdx.x<3)
	{
		if(threadIdx.x%2 == 1)
		{
			TriImgbk8[threadIdx.x] = TriImg9[threadIdx.x/2];
		}
		else
		{
			TriImgbk8[threadIdx.x] = (TriImg8[threadIdx.x] - TriMat8[threadIdx.x].z*TriImg9[threadIdx.x/2])/TriMat8[threadIdx.x].y;
		}
		
	}
	__syncthreads();
	
	if(threadIdx.x < 7)
	{
		if(threadIdx.x%2 == 1)
		{
			TriImgbk7[threadIdx.x] = TriImg8[threadIdx.x/2];
		}
		else
		{
			TriImgbk7[threadIdx.x] = (TriImg7[threadIdx.x] - TriMat7[threadIdx.x].z*TriImg8[threadIdx.x/2] - TriMat7[threadIdx.x].x*TriImg8[threadIdx.x/2 - 1])/TriMat8[threadIdx.x].y;
		}
	}
	__syncthreads();

	if(threadIdx.x < 15)
	{
		if(threadIdx.x%2 == 1)
		{
			TriImgbk6[threadIdx.x] = TriImg7[threadIdx.x/2];
		}
		else
		{
			TriImgbk6[threadIdx.x] = (TriImg6[threadIdx.x] - TriMat6[threadIdx.x].z*TriImg7[threadIdx.x/2] - TriMat6[threadIdx.x].x*TriImg7[threadIdx.x/2 - 1])/TriMat7[threadIdx.x].y;
		}
	}
	__syncthreads();

	if(threadIdx.x < 31)
	{
		if(threadIdx.x%2 == 1)
		{
			TriImgbk5[threadIdx.x] = TriImg6[threadIdx.x/2];
		}
		else
		{
			TriImgbk5[threadIdx.x] = (TriImg5[threadIdx.x] - TriMat5[threadIdx.x].z*TriImg6[threadIdx.x/2] - TriMat5[threadIdx.x].x*TriImg6[threadIdx.x/2 - 1])/TriMat6[threadIdx.x].y;
		}
	}
	__syncthreads();

	if(threadIdx.x < 63)
	{
		if(threadIdx.x%2 == 1)
		{
			TriImgbk4[threadIdx.x] = TriImg5[threadIdx.x/2];
		}
		else
		{
			TriImgbk4[threadIdx.x] = (TriImg4[threadIdx.x] - TriMat4[threadIdx.x].z*TriImg5[threadIdx.x/2] - TriMat4[threadIdx.x].x*TriImg5[threadIdx.x/2 - 1])/TriMat5[threadIdx.x].y;
		}
	}
	__syncthreads();

	if(threadIdx.x < 127)
	{
		if(threadIdx.x%2 == 1)
		{
			TriImgbk3[threadIdx.x] = TriImg4[threadIdx.x/2];
		}
		else
		{
			TriImgbk3[threadIdx.x] = (TriImg3[threadIdx.x] - TriMat3[threadIdx.x].z*TriImg4[threadIdx.x/2] - TriMat3[threadIdx.x].x*TriImg4[threadIdx.x/2 - 1])/TriMat4[threadIdx.x].y;
		}
	}
	__syncthreads();

	if(threadIdx.x < 255)
	{
		if(threadIdx.x%2 == 1)
		{
			TriImgbk2[threadIdx.x] = TriImg3[threadIdx.x/2];
		}
		else
		{
			TriImgbk2[threadIdx.x] = (TriImg2[threadIdx.x] - TriMat2[threadIdx.x].z*TriImg3[threadIdx.x/2] - TriMat2[threadIdx.x].x*TriImg3[threadIdx.x/2 - 1])/TriMat3[threadIdx.x].y;
		}
	}
	__syncthreads();
	
	if(threadIdx.x < 512)
	{
		if(threadIdx.x%2 == 1)
		{
			TriImgbk[threadIdx.x] = TriImg2[threadIdx.x/2];
		}
		else
		{
			TriImgbk[threadIdx.x] = (TriImg[threadIdx.x] - TriMat[threadIdx.x].z*TriImg2[threadIdx.x/2] - TriMat[threadIdx.x].x*TriImg2[threadIdx.x/2 - 1])/TriMat2[threadIdx.x].y;
		}
	}
	__syncthreads();

	pOutput[tiy*w+tix] += (unsigned char)TriImgbk[threadIdx.x];

}


__global__ void 
SobelTex( Pixel *pSobelOriginal, unsigned int Pitch, 
          int w, int h, float fScale )
{ 
    unsigned char *pSobel = 
      (unsigned char *) (((char *) pSobelOriginal)+blockIdx.x*Pitch);
    for ( int i = threadIdx.x; i < w; i += blockDim.x ) {
        unsigned char pix00 = tex2D( tex, (float) i-1, (float) blockIdx.x-1 );
        unsigned char pix01 = tex2D( tex, (float) i+0, (float) blockIdx.x-1 );
        unsigned char pix02 = tex2D( tex, (float) i+1, (float) blockIdx.x-1 );
        unsigned char pix10 = tex2D( tex, (float) i-1, (float) blockIdx.x+0 );
        unsigned char pix11 = tex2D( tex, (float) i+0, (float) blockIdx.x+0 );
        unsigned char pix12 = tex2D( tex, (float) i+1, (float) blockIdx.x+0 );
        unsigned char pix20 = tex2D( tex, (float) i-1, (float) blockIdx.x+1 );
        unsigned char pix21 = tex2D( tex, (float) i+0, (float) blockIdx.x+1 );
        unsigned char pix22 = tex2D( tex, (float) i+1, (float) blockIdx.x+1 );
        pSobel[i] = ComputeSobel(pix00, pix01, pix02, 
                                 pix10, pix11, pix12,
                                 pix20, pix21, pix22, fScale );
    }
}

extern "C" void setupTexture(int iw, int ih, Pixel *data, float *plogdata, int Bpp/*, float *foutput*/)
{
    cudaChannelFormatDesc desc;
    
    if (Bpp == 1) 
	{
        desc = cudaCreateChannelDesc<unsigned char>();
    } 
	else 
	{
        desc = cudaCreateChannelDesc<uchar4>();
    }

    cutilSafeCall(cudaMallocArray(&array, &desc, iw, ih));
    cutilSafeCall(cudaMemcpyToArray(array, 0, 0, data, Bpp*sizeof(Pixel)*iw*ih, cudaMemcpyHostToDevice));


	desc = cudaCreateChannelDesc<float>();
	cutilSafeCall(cudaMallocArray(&farray, &desc, iw, ih));
	cutilSafeCall(cudaMemcpyToArray(farray, 0, 0, plogdata, sizeof(float)*iw*ih, cudaMemcpyHostToDevice));
	
}

extern "C" void deleteTexture(void)
{
    cutilSafeCall(cudaFreeArray(array));
}


// Wrapper for the __global__ call that sets up the texture and threads
extern "C" void sobelFilter(Pixel *odata, int iw, int ih, enum SobelDisplayMode mode, float fScale, float *foutput) 
{
    cutilSafeCall(cudaBindTextureToArray(tex, array));
	cutilSafeCall(cudaBindTextureToArray(ftex, farray));

	
	//float *SNR;
	cutilSafeCall(cudaMalloc((void**)&SNR, sizeof(float)));
	cutilSafeCall(cudaMemset(SNR, 0, sizeof(float)));
	//CUDA_SAFE_CALL(cudaMalloc((void**) &d_a, sizeof(int)));
	

    switch ( mode ) {
		case  SIND:
		{
			
			dim3 snrthreads(BlockSize, BlockSize);
			SIND_SNR<<<1,snrthreads>>>(iw, iw, ih, 100, 100, 116, 116, SNR);
			cutilSafeCall(cudaMemcpy(foutput, SNR, sizeof(float), cudaMemcpyDeviceToHost));
		}
		break;
        case  SOBELDISPLAY_IMAGE: 			
            SobelCopyImage<<<ih, 384>>>(odata, iw, iw, ih, fScale );
            break;
        case SOBELDISPLAY_SOBELTEX:
            SobelTex<<<ih, 384>>>(odata, iw, iw, ih, fScale );
            break;
        case SOBELDISPLAY_SOBELSHARED:
        {
            dim3 threads(16,4);
#ifndef FIXED_BLOCKWIDTH
	          int BlockWidth = 80; // must be divisible by 16 for coalescing
#endif
        		dim3 blocks = dim3(iw/(4*BlockWidth)+(0!=iw%(4*BlockWidth)),
                               ih/threads.y+(0!=ih%threads.y));
        		int SharedPitch = ~0x3f&(4*(BlockWidth+2*Radius)+0x3f);
        		int sharedMem = SharedPitch*(threads.y+2*Radius);

        		// for the shared kernel, width must be divisible by 4
        		iw &= ~3;

        		SobelShared<<<blocks, threads, sharedMem>>>((uchar4 *) odata, 
                                                        iw, 
#ifndef FIXED_BLOCKWIDTH
                                                        BlockWidth, SharedPitch,
#endif
                                                		    iw, ih, fScale );
        }
        break;
    }
	
	cutilSafeCall(cudaFree(SNR));
	cutilSafeCall(cudaUnbindTexture(ftex));
    cutilSafeCall(cudaUnbindTexture(tex));
}
